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I. INTRODUCTION 

This paper deals with a new set of generalized Wannier functions (WF) which we 
argue has many advantages as a basis for large scale electronic structure calculations. The 
term Wannier functions refers to a complete and orthonormal set of local basis functions 
for representations of scalar wave functions. An example is the set of WF which span the 
same function space as the Bloch functions for a given energy band, or set of energy bands, 
in a crystal. Though they are often useful as a theoretical tool, energy band WF have 
seldom been practical in quantitative electronic structure calculations. There are several 
reasons for this: One difficulty is the numerical complexity of their construction; they are 
not the solutions of any simple differential equation and often must be constructed from non- 
orthogonal trial functions. Second, the advantage of their mutual orthogonality is offset by 
the existence of long range oscillatory tails which are necessary to ensure their orthogonality. 
Band WF are, however, exponentially localized provided there is a finite band-gap, and 
good localization is an essential requirement in quantitative calculations. One can define 
local Wannier functions which span the same function space as the plane waves in a given 
Brillouin zone, but such WF decay only algebraically and are of little practical use. 

"Generic" basis functions, i.e., functions with similar analytical structure, are partic- 
ularly advantageous in large scale electronic structure calculations. Examples of generic 
functions are plane waves and and gaussians, which are often used in electronic structure 
calculations. Evidently, the ease of computing Hamiltonian matrix elements with such an- 
alytical bases outweighs their disadvantages, namely, non-locality in the case of plane waves 
or non-orthogonality in the case of gaussians. Generic basis functions would appear to be 
especially desirable for vector or parallel computations. In contrast, an important disad- 
vantage of energy band WF is that they are not generic; their form depends on a given 
atomic environment and hence, they must be recomputed in the course of self-consistent 
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calculations. 

When the size of the basis becomes large, however, it becomes computationally advan- 
tageous to trade simplicity of matrix element construction in favor of a smaller, if more 
complicated basis. Sparseness of the Hamiltonian matrix is also a significant advantage. 
Yet another consideration is completeness. In the case of gaussians, it is not always clear 
whether the largest bases that can be implemented are sufficiently complete to cover the 
Hilbert space of interest. On the other hand, gaussian basis functions may be centered on 
arbitrarily positioned atoms, allow for spherical symmetry, and can represent well the cusps 
in atomic wavefunctions due to the divergence in the Coulomb potential. Mixed basis sets, 
e.g., a combination of gaussians and plane waves, have also been used. In general, such non- 
orthogonal bases can become overcomplete, as the finite basis is almost linearly dependent, 
and care must be taken to avoid numerical difficulties. 

These considerations suggest that desirable basis functions for large scale electronic struc- 
ture calculations should be orthogonal, complete, generic, local and possess useful analytic 
properties. To minimize basis set size, their locality should apply both to position and mo- 
mentum space. A set of functions which has all of these desired properties are the generalized 
Wannier functions wi „(x) introduced by one of us. Here / and n are momentum and site 
labels, respectively. In this paper we discuss the potential applicability of these Wannier 
functions, hereafter referred to as "phase space Wannier functions" (PWF), to electronic 
structure calculations. Like plane waves, the PWF form a complete set and have simple 
matrix elements with respect to the kinetic energy operator p = — V , thereby simplifying 
the calculation of Hamiltonian matrix elements. Like gaussians, they are well localized in 
position space, and matrix elements with respect to simple functions, e.g., 1, x, x , . . . may 
be stored. Moreover, they are shown below to be localized exponentially in both position 
and momentum space. This leads to a sparse Hamiltonian matrix and an overall basis set 
size which can be significantly smaller than that with plane waves. These properties suggest 
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that the PWF functions have significant advantages in large scale electronic structure cal- 
culations, and we argue that they are a potential replacement or supplement to plane wave 
and gaussian basis sets. These considerations have been checked by carrying out a detailed 
one-dimensional, model calculation with the PWF with convergence checks and comparisons 
with other approaches. However, a full scale test of the basis, e.g., with pseudopotentials 
and self-consistency remains to be done. 

The paper is organized as follows: In section II the basic properties of the functions 
wi j^{x) are given, together with an algorithm for their construction and a demonstration 
of their exponential localization. Some generalizations are also discussed. In section III 
the utility of these Wannier functions is illustrated with a simple one-dimensional model. 
Criteria for choosing the basis set size and lattice constant are determined, and convergence 
of the ground state energy with respect to basis size is evaluated. Section IV contains a 
summary and conclusions. 



II. PROPERTIES 

In this section we summarize the properties of the phase space Wannier functions (PWF) 
introduced in Ref. 1. In addition some generalizations are introduced and their localization 
properties in position and momentum space are established. 

A. Definitions 

By construction the PWF wi ^{x) are uniform translations of real functions wi{x), lo- 
calized about the sites of a regular lattice a; = na in position space 

wi^nix) =wi{x-2TTn). (2.1) 

Here a is an arbitrary lattice constant, which for convenience is set to 27r, / is a non-negative 
integer momentum index, as discussed below, and n is an integer site index; the functions 
wi{x) are of even(odd) parity for / even(odd). For other choices of a, normalized PWF are 
given by scaling the above functions: 

1 /9 

w;,n(a^;a) = (27r/a) / wi^nC^nx/a). (2.2) 

The functions wi „ are orthogonal and normalized such that their inner product is a delta 
function: 

{w^nlwi^n') = kvKn'- (2-3) 

The set of wi „(x) constructed in Ref. 1 are illustrated in Fig. 1(a) for < / < 2. 
In momentum space, Eq. (2.1) becomes 

*;,n(p) = e-^^^'^^^Kp). (2.4) 

where wi „(p) is the Fourier transform of wi n{x). The function wi{p) is peaked bimodally in 
momentum space, as illustrated in Fig. 1(b). More precisely, w;(p) is a real(pure imaginary) 

5 



function of even(odd) parity for even(odd) /, and is peaked near p = ±^(2/ + 1) (except 
for / = which is peaked only at the origin). Each function wi{p) roughly covers the 
disjoint intervals [^l, ^{l + 1)] and [—l{l + 1), —^l] in p, decaying exponentially rapidly into 
neighboring intervals. Thus one may imagine that phase space [— oo < x < oo, — oo < p < oo] 
may be spanned in terms of "fuzzy" rectangular blocks each represented by a function wi ^ 
either in position or momentum space (Fig. 2), whence our description "phase space Wannier 
functions." This localization property permits the basis to be tailored at a given site n simply 
by increasing the number of momentum states lmax{n) that are used without altering the 
form of the functions, changing the lattice constant a, or warping space. 

Several additional properties are needed to define the PWF uniquely. For example the 
functions of Ref. 1 are defined to have simple matrix elements with respect to the kinetic 
energy operator p = —d/dx, 

{"; (^n,n'-l - ^n,n'+l)> I = l' + l] 
Pl,n~n'^ I = ^'; (2.5) 

"r (<^n,n'+l - ^n,n'-l)) I = l' - I] 

where a; = //(^^r) and /3; „ decays exponentially rapidly with \n\. The /3; ^ are computed 
numerically in momentum space using Eq. (2.5). Choices of dispersion other than p will be 
discussed below. 

The three-dimensional generalizations of these functions are, as in the case of plane 
waves, Cartesian products, which are labelled by vector indices, / and n, 

Wf;.{f) = Wn^^l^{x)wn^^l^{y)wnj^{z). (2.6) 

The product form of Eq. (2.6) can considerably simplify the computation of matrix elements; 
for example, matrix elements of separable operators such as exp(— Ar ) factor into one di- 
mensional products. Thus matrix elements of the Coulomb potential can be computed as 
a single one-dimensional integral of a product of three one-dimensional matrix elements by 
use of the identity (1/t) = (27r^ ' ) Jq°° (i/T:exp(— /t r ). 
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The functions wi ^(x) have a number of other useful analytical properties. For example 
the coefficients q ^ in the Wannier function expansion, 

f{x) = ^ci^nWi^n{x), (2.7) 

l,n 

of a function f{x) with a Taylor expansion about x = na are related to the Taylor series 

coefficients of f{x) at na as follows: 

co,n =a^ (/M + ir(an)(|-)2 + 0{f(''^\na)a^)) 

cin =a^ (-/'(«^)(;f ) + 0{f"\an)a^) (2.8) 



'27r' 
. a ^ 



C2,n=a^(2r(an)(;^)2 + 0(/(-)a4)). 



Equation (2.8) follows from the relation between the moments of wi jj_{x) and the derivatives 
of wi n{p) at p = 0, as discussed in Appendix A. 

As a consequence, the expansion of the constant function requires only the / = PWF 
(at all sites); the expansion of a linear function requires only the / = 0, and 1 functions (at 
all sites); and so on. This result holds locally as well: in the neighborhood of a site x = na 
the first two PWF's (/ = 0, 1) adequately represent an arbitrary function as long as it does 
not vary rapidly on the scale of a. For sufficiently narrow cells, these functions therefore 
serve as a complete set of "shape functions," as in the finite element method of Ref. 5. While 
the shape functions of the finite element method can be more compact, their finite range 
leads to poor convergence in momentum space; also, these shape functions are polynomial 
complete only to some low order. 

B. Algorithm for Construction of wi ^ 

An algorithm together with a Gibbs description of a computer code for constructing 
PWF wi ^ with the above properties has been given in Ref. 1. In this Section we present an 
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alternative procedure which has the following advantages: The Wannier functions generated 
can be shown to form a complete and orthogonal basis, and the new procedure can be 
generalized to make use of dispersions other than p . Lastly, a recursion relation for the 
PWF is discussed. 

To prove completeness, and orthogonality of the PWF basis, we show that the PWF's 
are the eigenfunctions of a Hermitian matrix H representing the operator p = —d/dx. 
The diagonal matrix elements of H, which are determined from the Pi „, are not known 
beforehand, but are computed by an iterative method. 

The site-translation property of Wannier functions allows one to work in a reciprocal 
space description in which p = m + q with m an integer, and q in the first Brillouin zone. In 

this picture the orthogonality relation, Eq. (2.3) becomes 

+00 

E(^^,0kl>)e^"''^" = E w^{rn + q)wiim + q) = Su>. (2.9) 

n m——oo 

Similarly, defining Hiii{q) = Sn(w; gb l""^!' n)c *'^"') the p matrix elements, Eq. (2.5), 

become 

+00 ( 2iai sin27rg, / = /' + !; 

Hi,l'iq)= Yl wlirn + q)wi,im + q)im + qf=!pi{q), 1 = 1'; (2.10) 

m=-oo [ —2iaii sin27rg, / = /' — !; 

with ai = //(47r), and A(g) = S„/3;,^e2«^9. 

For each g, the tridiagonal matrix Hi i/{q) defines an eigenvalue problem for which the 

eigenvalues are known explicitly, and the eigenfunctions are the desired PWF. 

00 

X] Hi^i'{q)wi'{m + q) = {m + qfwi{m + q). (2.11) 

l'=Q 

The diagonal elements of Hi ii{q) are computed numerically as follows: First, the range of / is 
made finite by limiting it to {0, ■ ■ ■ ,L} with L even. The set of eigenvalues is also truncated 
to {i~^^ + ?) 5 ("5-^ + 1 + q) , ■ ■ ■ , i^L + q) } to reflect the truncated momentum space. 
The truncation introduces signiflcant error in quantities near the upper end of momentum 



space; for example, the values of Pii^q), wiim + g), and W{){\L + q) computed in this 
approximation differ significantly from their true values. On the other hand, the effects of 
truncation diminish rapidly as one moves away from the upper end, and accurate Pi{q) can 
be computed for / small compared to L. 

We evaluate the characteristic equation for each of the L + 1 eigenvalues to obtain a set 
of equations determining the L + 1 unknowns Pi{q), 

det\Hi^l.{q) - (m + g)2<5^^H = 0, for m = {-IL, ■ ■ ■ , iL}. (2.12) 

These equations may be solved by a fixed-point iterative scheme as each Pi{q) may be 
readily expressed in terms of the others because the determinant is linear in each (3i{q). The 
expected bimodality of the wi{p) suggests that the starting values for the Pi{q) be chosen to 
be near (i/ + q) for / even, and {—^{l + 1) + q) for / odd. Once the /3;(g) are known, a 
quick diagonalization yields the w(m + g). 

Note that the method of Ref. 1 was to solve for the wi{p) without computing the /3;(g) 
first. This is more efficient and is the preferred method for generating the PWF basis; 
however, the present prescription can also be generalized to dispersions other than p as 
discussed in Section II C, below. 

Recursion relations between the functions wi{p) are also useful. From the matrix elements 
in Eq. (2.5) one has 

ai+iwi+iip) = aiwi_i{p) + i^-lAM*Kp), (^ > 0, ao = 0), (2.13) 

li sm iTip 

Thus, given the function wq{p), the remaining functions and the Pi{p) can be determined 
in principle by recursion. Note that the functions wi{p) have zeros at the half integers (at 
which the denominator 2isin27rp vanishes) with the exception oip = ±5/, and p = ±i(/ + l) 
(at which the numerator p — (3i{p) also vanishes). Tables of the functions wi{p) are available 
from the authors. 
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C Generalization to Arbitrary Dispersion 

An advantage of the construction algorithm outlined in Section 11 B is that the general- 
ization to other dispersions H = e{p) is easily carried out. We consider only dispersions that 
do not have the periodicity of the reciprocal lattice and are even functions of p. Under these 
restrictions the generalized WF can be chosen to have definite parity and cover the phase 
space in the same manner as before. 

Matrix elements of the dispersion are again required to have a tridiagonal form: Eqs. 
(2.10) and (2.11) are replaced by 

+00 ( 2iai sm27cq, / = /' + 1; 

Hll'iq) = Yl *n^ + ?)*;'(^ + Q) <rn + q) = < /3i{q), I = l'; (2.14) 

m=-oo [ —2iaii sin27rg, / = /' — 1; 

with /3;(g) = Sn/3;,ne2'^*"'^ and 

Yl Hi^VWi,{m + q) = e{m + q)wi{m + q). (2.15) 

V=Q 

By an analysis similar to that of Ref. 1 we find that the coefficients a; are given by e'(i/)/47r, 

where e'(p) = de{p)/dp. The remaining steps of the construction procedure are carried out 

as in the previous section. The iteration is started with the /3/(g) chosen to be near e(i/ + q) 

for / even, and e( — i(/ + 1) + g) for / odd. The effect of dispersion on the PWF is illustrated 

in Fig. 3, where functions for the dispersions e{p) = p and e{p) = p are compared. 

D. Exponential Localization 

The Wannier functions wi ^j of Ref. 1 appear to be localized exponentially in position 
space and their fourier transforms in momentum space. This is surprising in view of a 
theorem proved independently by Balian and by Low, which shows that it is impossible to 
construct Wannier functions which are exponentially localized about the sites of a regular 
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lattice in phase space with cell area 2tt. In this section we discuss how this theorem can be 
evaded to construct functions which do possess exponential locahzation. 

The non-existence theorem of Ref. 8 applies specifically to functions which are uniform 
translations in phase space of a singly peaked function f{p), as in the functions of von 
Neumann, 

^;,n(p) = /(p-0e-^2"^", (2.16) 

for integers l,n. As a consequence expansions in "coherent states," where f{p) is a gaussian, 
converge very slowly {i.e., algebraically), so their usefulness in quantitative calculations is 
limited. By contrast, the PWF are bimodal in phase space. They are composed of functions 
peaked about p = ±^(2/ + 1), as in 

9l,n{P) = [f(P - i(2^ + 1)) + (-I)V(-P - i(2/ + 1))] e-2.^^ (2.17) 

for integers / > and n. This bimodal form is sufficient to render the theorem in Ref. 8 
inapplicable. Strictly speaking, the bimodal form of gi j^{p) implies that the rms width of 
such functions in momentum space is 0{l) for large /; but the two peaks f{p) in Eq. (2.17) 
are each well localized about |p| = ^(2/ + 1). 

To see how bimodality of gi „ leads to exponential localization, we consider the problem 
of generating a complete set of orthonormal, Wannier functions wi „ starting from trial 
functions of the form of Eq. (2.17), where f{p) is a normalizable, even function assumed 
to be well localized about the origin. The functions wi ^ are constructed by a symmetric 
orthogonalization of the functions gi ^, using the inverse square root of the overlap matrix 

Sln,l'n' = {9l,n\9V,n') '■ 



/jS I )ln,l'n'9l',n'- (2.18) 



l',n' 

The matrix S~ ' can be calculated by going to a diagonal representation in terms of the 
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solutions to the eigenvalue problem, 

Z2 ^ln,l'n' Cl',n' = -^Q,n- (2-19) 

l',n' 

This construction is carried out in Appendix B where it is shown that ("S*^ / )^^ ^/jj/ decays 
exponentially both in \n — n'\ and in \l — l'\ for large /. This result implies exponential 
localization with the same decay constants for the Wannier functions, wi ^ defined by 
Eq. (2.18), for sufficiently localized trial functions f{p) {e.g., gaussians). Moreover, the 
translational invariance (iS*^ ' );^ ^/^/ in / for large / implies that the Wannier functions in 
momentum space will be nearly identical for large /, a form which is already evident from 
Fig. 1(b). 
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III. APPLICATION TO QUANTUM 
MECHANICS IN ONE-DIMENSION 



A. Schrodinger Equation in the PWF Basis 

As an illustration of their properties and the application of the PWF we use them as a 
basis for solving the one-dimensional Schrodinger equation 



Hip 



+ Vix) 



iIj = EiIj. (3.1^ 



dx'^ 
We take V{x) to be an attractive single- well potential of the form 

V{x) = -A sech2(a;), (3.2) 

with A = 35/4 (Fig. 4). For this choice, the bound state solutions iIj{x) are known analyt- 
ically, so that precise comparisons are possible: there are three bound states of which the 
lowest is at Eq = -25/4, and ipQ = (8/37r)V2sech^/2(a;). 

In the PWF basis, ipi^x) = T^i^ci^wi^^x), and the Schrodinger equation (3.1) is equiva- 
lent to the matrix eigenvalue problem, 

/ , Hin^ii^i Cyni = E Ci^. (3.3) 

I'n' 

It is convenient to compute the kinetic energy part of H in momentum space using Eq. (2.5) 
and the potential energy part in position space, so that 

The sum over n' in Eq. (3.3) is restricted to box of length L so that —^L<x = na< \L, 
for a given choice of L and lattice constant a, while that over /' is restricted to /' < lmax{i^')- 
The exact coefficients q „ can be obtained conveniently in momentum space, 

ci,n = / dpwl^{p)ilJo{p), (3.5) 
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where V'o(p) = (8/37r)(2/37r) / |r(| + ^ip)\. These coefficients are given in Table 1 for 
lattice constants a = 0.5, 1.5, and 3.5, and a fixed overall length L = 10.5, a value chosen 
to be a few times the width of the potential well V{x). Choosing, for example, 0.01 as 
a cutoff for the phase-space cell occupations |c;„| , the representation of the ground state 
wave function requires 5 PWF, i.e., (/ = 0; n = 0, ±1,±2), for lattice constant a = 0.5; 3 
PWF, i.e., (/ = 0;n = 0,±l), for a = 1.5; and just 2 PWF, (/ = 0,2,n = 0), for a=3.5. For 
comparison, a plane wave basis to the same 1% accuracy requires 9 plane waves, k = 27rm/L 
with —4 < m < 4 (see Table 2). With the cutoff set at 0.003, the corresponding numbers 
are 7 Wannier functions for a = 0.5, 5 for a = 1.5, 6 for a = 3.5, and 11 plane waves. 
Note that the example we have chosen is centered at x = 0, so some coefficients vanish 
automatically; similarly, if symmetrized plane waves were used, only 5 (or 6 for a cutoff of 
0.003) even functions would be needed. However, in practice one cannot assume symmetry. 
Without symmetry, the n = 0, / = 1 amplitudes need not vanish, increasing the counts for 
the Wannier functions by one for all three values of a; similarly one needs all plane waves. 
Note also that the number of plane waves needed depends on the length of the system L, 
while the number of Wannier functions is determined locally and is independent of L. In 
the present example, the number of plane waves can be reduced by cutting down L, but this 
restricts the location of the orbital. 

B. Phase Space Criteria for Basis Set Size 

We now discuss the dependence of the coefficients q „ on / and n from phase-space 
considerations. To the extent that a semi-classical description is valid, these coefficients 
should be small for all phase space cells which do not overlap the classical trajectory. This 
indeed appears to be the case and suggests a criterion for adopting a given lattice constant 
and basis set size. One choice of the lattice constant a is that value for which the "coverage" 
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"^000 \^l n\ is optimal, i.e., the value for which the minimum number of phase space cells 
overlap the classical region of occupied states. From Table 1 the coverages at a cutoff of 0.01 
in IQj^I are 0.978, 0.985 and 0.975 for lattice constants a = 0.5, 1.5, and 3.5, respectively. 
For a cutoff of 0.003, the corresponding figures are 0.990, 0.999 and 0.997. For both cutoff 
choices, the value a = 1.5 maximizes the coverage with a minimum number of functions and 
is therefore optimal for the values being considered. Inspection of Fig. 5 indicates that this 
lattice constant corresponds roughly to the size of the phase space orbit. Were the lattice 
constant chosen to differ greatly from this value, the phase space cells would not "tile" the 
orbit efficiently, and more PWF's would be needed to achieve comparable coverage of the 
ground state wavefunction. 

For comparison, a plane wave expansion spans phase space with elongated cells of width 
Ax = L and height Ap = 2tt/L, and therefore requires 2prnaxL/2'K states, for an equivalent 
description. The value of pmax is related to the depth of a potential well. For example, 
in the case of a point charge Ze, ip{x) ~ exp(— Zx/co) and Pmax ~ ^'/oo, where Oo is the 
Bohr radius. The use of pseudopotentials, therefore, can considerably reduce the size of 
basis required for a given calculation, either in a plane-wave or a PWF basis. Moreover, one 
expects from general phase space volume considerations {i.e., from the uncertainty principle) 
that the number of Wannier functions needed will be of order 2j dxp{x)/27ih which is just 
the area of phase space within the classical orbit at energy E, in units of 2'n-h. The Wannier 
functions are therefore expected to be more efficient than the plane waves by a factor roughly 
equal to the ratio of these phase-space areas (phase-space volumes in three dimensions). It 
is also of interest to examine the convergence of the "coverage" vs. total basis size. For 
smooth potentials, the convergence of the coverage with respect to / at a given n appears 
to be exponentially rapid, which is what one expects based on WKB arguments. For the 
one-dimensional example discussed above, it is observed that the number of plane waves 
required to give an accuracy of 1% (or 0.3%) in the coverage is, ignoring symmetry, about 
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9/6 (or 11/6) times that for the Wannier functions; in three dimensions the corresponding 
factor would be cubed, i.e., with the plane wave basis larger by a factor of 3 to 6, depending 
on the choice of cutoff. Thus because it can be adjusted locally, the PWF basis can be 
significantly smaller than a plane wave basis. 

C. Convergence of Ground State Energy 

Finally, we have computed the coefficients q ^j and the eigenvalue Eq for the lowest state 
by straightforward diagonalization of the Hamiltonian matrix in Eq. (3.3) for several different 
basis set sizes. The results are shown in Table 3. Note that convergence of the ground state 
energy and the convergence of the coverage are comparable since both vary quadratically with 
wavefunction error. For example, for the case a = 1.5 with 5 basis functions the difference 
from unity of the coverage, i.e., 1 — S|q ^\ , is 0.0027, and the relative error in the ground 
state energy, AEq/Eq, is 0.0008. For a = 3.5 with 6 basis functions, the corresponding 
numbers are 0.0046 and 0.0102, respectively. This also indicates the importance of total 
coverage as a gauge of convergence of a given basis set. 
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IV. CONCLUSIONS 

We have suggested that phase space Wannier functions may be an attractive basis for 
large scale electronic structure calculations, as they constitute an orthonormal, complete, 
local and generic basis set. Their localization in position and momentum space permits the 
basis set to be optimized locally, and they have many convenient properties, e.g., simple 
short-range matrix elements of the kinetic energy operator. On the other hand, we have not 
constructed an analytical representation of these functions; at present, they are tabulated 
numerically, and typically 100 points are needed to represent their fine structure. Detailed 
properties of these functions and some generalizations have been given, together with an 
algorithm for their construction and proof of their exponential localization properties. 

An illustrative application to the Schrodinger equation in one-dimension has also been 
given. The results of our one-dimensional example suggest criteria for choosing the basis set 
size and lattice constant a based on semi-classical, phase space considerations. Compared 
with a plane wave basis set, the phase space Wannier functions permit a significant reduction 
in basis set size, since the basis size can be adjusted locally. Even a small reduction is 
significant in three dimensions, since the number of operations required in a straightforward 
matrix diagonalization varies as the cube of the basis size. 

The next step would be to carry out a three-dimensional electronic structure calculation 
using these functions, i.e., a calculation analogous to a comparable plane- wave, pseudopo- 
tential calculation. Important developments needed for self-consistent calculations are a) 
a good representation of pseudopotential matrix elements, and b) an efficient scheme for 
solving Poisson's equation in the PWF basis. 
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APPENDIX A 



Taylor and Wannier Function Expansion Coefficients 

In this appendix we discuss the relation between the Wannier function expansion coeffi- 
cients q „ of a function f{x), 

f(.)^j:c„..,,„M. e,,„^ /...,,,.;.)/(.). (Al) 

l,n 

and the Taylor series expansion coefficients of f{x) about the point x = na , i.e., 

fi.^) = E , (^ - n^V- (A2) 

i 

Using equations (Al), (A2) and (2.2), two expressions for the expansion coefficients are 

readily obtained. 

q,„ = (a/27r)^ Y. ^-^^^{a/2ny Jdxx^wi{x) 
J 

J 
= a^^Z!!lM(,,/2.).-p)(o) (A3) 

j 

where /i;j is j moment of wi{x), and wY {p) is the j derivative of wi{p). This last equa- 
tion follows from the equality between the moments and the derivatives of wi{p) evaluated 
at the origin. 

^^l^^ = {27,)U^w\'\Q). (A4) 

Now: the sum in Eq. (A3) simplifies considerably since (a) many of the moments vanish 
by symmetry arguments, and (b) many more vanish by the structure of the wi{p) at the 
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origin. From the recursion relation Eq. (2.13) one can determine the location as well as the 
order of the zeros of wi{p). In particular, at the origin the following derivatives vanish, 

^{•^)(0) = forj</. (A5) 

Further, several of the non-vanishing derivatives can be deduced from the recursion and the 
orthogonality equations: wq{0) = 1, Wq(0) = —1, u'j^(O) = i, wi^iO) = —4; the remaining 
derivatives, if needed, can be calculated numerically. Using Eq. (A3) and known values of 
the derivatives, we obtain the results given in Eq.(2.8). 
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APPENDIX B 



Exponential Decay of Phase Space Wannier Functions 

In this section we construct the inverse square root {S^ ' )ij^ yj^t of the overlap matrix 
^Inl'n' = {91 nidi' n') ^^T a set of trial Wannier functions gi ^^ given by Eq. (2.17). The 
computation is most conveniently carried out in a diagonal representation, i.e., in terms of 
the eigenstates q ^ of the matrix S'|„ |/„/: 

Z2 ^ln,l'n'Cl',n' = -^Q,n- (Bl) 

l',n' 

This calculation is carried out in several steps as follows. First by translational invariance 
in position space, we may take q ^ to be of Bloch form, c/ ^ = ci{k) exp{ikn) so that Eq. 
(Bl) reduces to a one-dimensional eigenvalue problem, J^iiSi ii{k)cir{k) = Xj^ci{k), where 
'S'n'(^) — ^nexp(27riA;n)S'^Q ^/jj. Second, note that by construction the form of the matrix 
Si ii{k) depends on whether / and /' are even or odd. Thus we define / = 2/i + a, (cr = 0, 1). 
We find below that for large /,/' the matrix Si ii{k) becomes translationally invariant in fi. 
The eigenvalue problem is thus analogous to a tight-binding problem on a semi-infinite chain 
with two basis states {a = 0, 1) per unit cell. The limiting form of Si ii{k) is seen to be 

Siv{k) ^ s,,,{k,ii- /) = [s^Ak.i^- /) + (-i)"+"'^,,,(-A;,/^ - /i')] , (B2) 

where nontranslationally invariant terms which vary as S^j fji{k, ji + jj!) have been neglected; 
those terms are negligible for sufficiently localized functions /. In terms of the functions 
f{p) in Eq. (2.17) the quantity in Eq. (B2) is 

~Sa,a'{k,li) = Y,f\k + m-\{a + l))f{k + m + iJi-\{a+l)). (B3) 

m 

Next, the eigenvalues of Eq. (Bl) can be obtained from the large / limit where, by transla- 
tional invariance in /x, ci{k) — ?■ exp{±i2TTqfi)ca{k, q) . Then one obtains the 2x2 eigenvalue 
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problem, 

Sa,a'{k, q)c^{k, q) = Xk,qCa{k, q), (B4) 

where 

Sa,a'{Kq) = J]5,^'(A;,/x)e^'^^ = S^^^^{k,q) + {-lf+-'s^^^,{-k,q). (B5) 

It is easily verified that the functions S'crcr' factor as S(jfji{k,q) = 9'^{k,q)9^i{k,q), where 

Oa{k, q) = J2 e'^'^^'^fik + m - i(2a + 1)). (B6) 

m 

The functions in Eq. (B6) are analogous to those discussed by Balian and by Low. They show 
that 6cr{k,q) must have at least one real zero in the Brillouin zone [0<A;<l,0<g<l]. 
This property is sufficient to rule out exponential decay for Wannier functions constructed 
from the singly-peaked trial functions of Eq. (2.16). For the functions 6(j{k,q) one can 
also establish the identity 9i{k,q) = exp(i27rg)^Q(— /c, g), from which it follows that the 
2x2 matrix S^^'{k, q) is diagonal, with diagonal elements given by \OQ{k, q)\ + \9i{k, q)\ . 
Hence, barring the unlikely possibility of common zeroes both in 6q and 6i, the eigenvalues 
are positive definite for real k, q. For gaussian trial functions f{p), for example, there is only 

1 "^ 

a single zero of 9Q{k, g) at g = ^ and k = ^. 

Finally, the inverse square root matrix is obtained by quadrature, and we obtain 

{S-'/\nl'n'= II dkdq V cf (fc g)c^(A:, g) ^^,^,(^_„,) 

''"''" JJBZ \^i[|^o(A:,g)|2 + |^l(A;,g)|2]l/2 

where ci{k,q) = Ai{k,q) cos{qfi + 6i{k,q)) are the exact eigenstates, 6i{k,q) being a phase 

shift which is / dependent for small /. Since the integrand has no zero for real k, q in the 

Brillouin zone, ("S*^ ' );^ ;/^/ will decay exponentially in \n — n'\ and in |/ — /'| for large /, at 

rates determined by the nearest zeroes of A (A;, g) in complex k, q space. Applying this result 

to Eq. (2.18), one sees that the Wannier functions w; „ have the same decay rates, i.e., 

there is decay constant hn such that exp{hn)wi ^ — )■ 0,72 — )• oo,h < hn, and similarly for 

decay with respect to /. In Ref. 1 these decay constants are found empirically to be hn = 2.9 

and h-i = 0.46, respectively. 
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FIGURE CAPTIONS 

Fig. 1. Phase space Wannier functions wi for the dispersion relation p with lattice constant 
a = 2tt (a) in position space: / = 0, solid line; / = 1, dotted line; and / = 2, solid-dashed 
line; and (b) in momentum space, plotted with respect to the central momentum (2/ + l)/4. 
For odd /, the imaginary part is plotted. Note that the functions are indistinguishable for 
/ > 2. 

Fig. 2. Schematic plot of phase space cells (disjoint for / > 0) corresponding to peaks in wi ^ 
for (a) even / and (b) odd /. Also shown at is the approximate shape of wi{p) in momentum 
space. 

Fig. 3. Effect of dispersion relation e{p) on the PWF; shown are functions wi{p) in momentum 
space for dispersion relation p (solid line) and p (dashed line): / = (upper curves) and 
the imaginary part of Wi{p) for / = 1 (lower curves), all with lattice constant a = 27r. 

Fig. 4. Potential V{x) = — Asech (x) with A = 35/4. Horizontal lines correspond to the 
three bound-state energy levels. 

Fig. 5. Classical phase space orbits p{x) = \/ E^ — V{x) of bound states of the potential 
Vix) = — Asech (x) with A = 35/4. Superposed are phase space cells for different choices 
of lattice constants: a) a = 0.5, b) a = 1.5, and c) a = 3.5. Tables of coefficients for these 
choices of a are given in Table 1. 
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TABLE CAPTIONS 

Table 1. Contribution to the ground state wave function |c/„| from the phase space cell 
l,n, for different choices of lattice constant: a) a = 0.5; b) a = 1.5; and c) a = 3.5. 

Table 2. Contribution to the ground state wave function |cy^| for a plane wave expansion in 
the interval -L/2 <x< L/2, with L = 10.5 and k = 27rm/L, for m = 0, ±1, ±2, ■ ■ ■ ± 6. 

Table 3. Ground state energy Eq for various basis set sizes {l,n}. 
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Table 1. 





|c/,nP 


l\n 





1 


2 


3 


a = 


0.5 





0.41400 


0.23226 


0.04984 


6 X 10-3 


1 


0.00000 


3 X IQ-^ 


1 X 10-3 


2 X 10-4 


2 


2 X 10-4 


3 X 10-5 


2 X 10-5 


2 X 10-5 


a = 


1.5 





0.92235 


0.03112 


1 X 10-6 


5 X 10-^ 


1 


0.00000 


0.00731 


4 X 10-5 


1 X 10-*^ 


2 


1 X 10^4 


1 X lO-'^ 


1 X 10-*^ 


1 X 10-*^ 


a = 


3.5 





0.84512 


7 X 10-4 


1 X 10-5 


3 X 10-8 


1 


0.00000 


5 X 10-3 


1 X 10-5 


2 X 10-8 


2 


0.12984 


6 X 10-3 


3 X 10-5 


7 X 10-8 



Table 2. 



A; 


Ck' 


0.0000 


0.24702 


0.5984 


0.20007 


1.1967 


0.11042 


1.7952 


0.04524 


2.3936 


0.01500 


2.9920 


0.00429 


3.5904 


0.00111 



Table 3. 





N 


il;n) 


Eo{l,n} 


a = 0.5 


5 


(0;0,±1,±2) 


-4.8479 


9 


(0;0,±1,±2) (1;±1,±2) 


-6.0017 


a= 1.5 


3 


(0;0,±1) 


-6.0119 


5 


(0;0,±1) (1;±1) 


-6.2333 


a = 3.5 


2 


(0;0) (2;0) 


-5.9969 


6 


(0;0) (1;±1) (2;0,±1) 


-6.1860 
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